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Recent proposals of measuring bipartite Renyi entropy experimentally involve techniques that 
hold exactly for non-interacting quantum particles. Here we consider the difference between such 
measurements and the actual Renyi entropy for ground state fermion ab initio molecular systems. 
To calculate various entanglement measures we extended several different techniques for use with 
variational Monte Carlo, Hartree-Fock and quantum chemistry methods. Our results show that 
in systems with strong electron correlations the Renyi entropy may not be accurately determined 
with the proposed measurements. In addition, we find significant physical insight can be gained 
by calculating entanglement properties in molecular systems. We see that the Renyi entropy and 
the entanglement Hamiltonian encode information about the character of the covalent bonds in a 
molecular system and that such information may lead to better descriptions of bonded systems that 
have been traditionally hard to describe with standard techniques. In particular our results for the 
C2 dimer suggests all eight valence electrons play a significant role in covalent bonding. 



The study of bipartite entanglement entropy and 
other entanglement properties has led to new in- 
sights in electronic properties of extended systems 
as seen in many recent studies [1-11]. Real space 
bipartite entanglement properties are derived from 
reduced density matrices in which space is divided 
into two regions and spatial degrees of freedom are 
traced out over one of these regions. This suggests 
they contain information of how quantum systems 
are spatially connected to each other. Interpreting 
entanglement properties to develop an understand- 
ing of how they are useful in describing quantum 
systems is still an active field. This is in part due 
to the relatively small number of studies of entan- 
glement properties in interacting systems that are 
not related to fractional quantum Hall studies or ID 
Hamiltonians[12-15]. Understanding entanglement 
properties of interacting systems is of prime interest 
in recent proposals that suggest direct measurements 
of entanglement and Renyi entropies are experimen- 
tally possible[16-20]. 

We consider two specific experimental proposals 
for measuring entanglement properties in this work: 
the correlation function approach[21-23] and parti- 
cle number fluctuations[16, 17]. It has been shown 
that both of these methods are exact when the 
ground state wave function of interest represents a 
non-interacting system. However it is unknown what 
biases these proposed measurements will have in de- 
termining entanglement for interacting systems. On 
a more fundamental level we consider the relation- 
ship between entanglement properties and the char- 
acterization of molecular systems. In particular we 
explore the idea that the reduced density matrix, 




FIG. 1. The bipartite entanglement properties are cal- 
culated on space partitioned into two half spaces for 
homonuclear diatomic molecules. Region A can be con- 
sidered the left half space, and region B, which is traced 
out, is the shaded region on the right. The correlation 
method uses a large grid in region A, while the swap 
and the fluctuation methods are sampled continuously 
in region A to compute the Renyi entropies. 



generated by tracing out real space degrees of free- 
dom, contains information about the bonding prop- 
erties of a system. In this work we develop the pre- 
liminary framework for using these tools to analyze 
covalent bonds. 

The Renyi entropies are among several quantities 
of interest that are related to the reduced density 
matrix of a subsystem. The reduced density matrix 
that describes a subsystem will in general have finite 



Renyi entropies which are given by 

S„(A) = -J-]np&((p A )»)] (1) 
1 — n 

where n > and p A refers to the reduced density 
matrix calculated when the degrees of freedom out- 
side of region A are traced out from the density ma- 
trix p of the entire system. The Renyi entropies 
for integer values n > 1 are directly accessible from 
quantum Monte Carlo (QMC). In this work we are 
primarily interested in S2 and all further references 
to Renyi entropy will refer to S2 except when explic- 
itly stated otherwise. 

We use three techniques for calculating the Renyi 
entropy, which we refer to as the swap method, cor- 
relation method and fluctuation method. The swap 
method is an unbiased estimator, and is the bench- 
mark for which the other two methods should be 
compared. The other two methods involve calcu- 
lating intermediate values that are closely related 
to quantities that in principle can be measured 
experimentally. Through a combination of vari- 
ational Monte Carlo (VMC) and quantum chem- 
istry methods we can compute the Renyi entropy 
of the fluctuation and correlation methods and com- 
pare to the value that can be computed with the 
swap method[12]. The quantities to be calculated 
in the particle fluctuation method are the moments, 
/i TO , of the particle number operator in region A, 
Mm = (*|(iV A ) m )|*)/(*|*). With a finite number 
of moments it is possible to estimate S n [16]. Calcu- 
lating these moments can be performed directly in 
VMC. 

The correlation method is performed in a quite 
different manner and involves the expectation value 
of specific single body operators. The expecta- 
tion value of these operators can be calculated 
with the one particle reduced density matrix 71 us- 
ing Hartree-Fock and configuration interaction tech- 
niques. An important quantity that is calculated 
as an intermediate step in the correlation method is 
the entanglement Hamiltonian, Hent, which is de- 
fined as pa — ree _Wen( . The normalization k is set 
to ensure Ty(pa) = 1. All three methods are exact 
for Hartree-Fock wave functions. The Hartree-Fock 
results are shown in figure 3 for the C2 dimcr. The 
details of these methods are described in the supple- 
mental material. 

The systems considered in this work are those 
of the homonuclear first row dimers with ground 
state wave functions that are known to be spin 
singlets [24]. For these Hamiltonians we consider a 
partition of two half spaces as shown in figure 1 . Re- 
gion A consists of a half space containing one atom, 
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FIG. 2. The bipartite Renyi entropy, Sa, calculated for 
several diatomic molecules at their equilibrium geome- 
tries with the three different techniques described in the 
text. The swap operator is an unbiased estimate while 
the fluctuation and correlation methods are what would 
be measured with experiment for the same ground state 
wave functions. 

and region B consists of the half space containing the 
other atom. The Renyi entropy of several dimers at 
their equilibrium geometries is shown in figure 2 us- 
ing the three methods. In figure 3 we plot the Renyi 
entropy as a function of atomic separation for the 
C2 dimer over a range in which the covalent bonds 
between the atoms are stretched and broken. 
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FIG. 3. S2 Renyi entropy calculated for the C2 dimer for 
the multi-determinant wave functions. The Hartree-Fock 
values are plotted as well for reference. The equilibrium 
distance is at 1.2 A. 

For these systems, Hartree-Fock wave functions 
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are not always accurate representations of the 
ground state wave function. We have tested a vari- 
ety of different wave functions and determined that 
multi-determinant wave functions using a range of 
100-1000 determinants, depending on the system, to 
be quite accurate [24-26]. These wave functions are 
sufficiently compact to include all determinants in 
the VMC calculation, and in this form it is straight- 
forward to compute all three techniques for calculat- 
ing the Renyi entropies. 

The results in figures 2 and 3 show that the corre- 
lation and fluctuation methods capture some of the 
qualitative behavior of the Renyi entropy but not 
in the limit of large atomic separation of C2. It is 
interesting that these two methods consistently un- 
derestimates or overestimates the Renyi entropy for 
all the systems considered here. Neither estimate is 
consistently accurate; however when considered to- 
gether they appear to bracket the actual value of the 
Renyi entropy. This can be partly understood from 
known properties of the fluctuation method. It has 
been shown that estimates of Si, from the fluctu- 
ation method are strictly less than or equal to the 
actual value of Si [19]. We observe this bracketing 
holds for all Renyi entropy calculations considered 
in this work, and it would certainly be of interest to 
explore its generality further. 

In the region of the equilibrium geometry for C2 
the Renyi entropy of the multi-determinant wave 
function is in relatively good agreement with the 
Renyi entropy of the Hartree-Fock wave function. 
These results show the Renyi entropy increases as 
the atoms are squeezed past the equilibrium ge- 
ometry followed by a maximum and reduction in 
the entropy. The agreement between the multi- 
determinant and the Hartree-Fock wave functions 
breaks down as the separation between the atoms 
increases, which is also the region in which the fluc- 
tuation and correlation methods show large devia- 
tions from the benchmark swap results. The Renyi 
entropy from the fluctuation method decays quickly 
when compared to the swap operator. In contrast, 
the correlation method does not allow the electrons 
to disentangle in the bond breaking region. In this 
limit it is expected that the fluctuation method sig- 
nificantly underestimates the Renyi entropy as elec- 
tron number fluctuations between subsystems are 
only part of the Renyi entropy for an interacting 
system. Two regions can have finite entanglement 
without fermions exchanging between them such as 
in a dispersion bonded system. 

Interpreting the value of the Renyi entropy can 
lead to insight into properties of molecular systems 
and here we consider how it might be related to the 



bonding properties. We consider the cases of C2 
and N2 at their equilibrium geometries, where the 
Renyi entropy is slightly greater than 4.0 for both 
molecules, and we ask what is the relevance of the 
magnitude of these Renyi entropies. We develop an 
interpretation of these values by considering the en- 
tanglement of a single non-interacting electron di- 
vided evenly in half by region A and B for example 
the case of a particle in a box that is split by the 
midpoint of the box. The results of this model is a 
density matrix p^ that has two eigenvalues (1/2,1/2) 
which yields a value ln(2) for the S2 Renyi entropy. 
For our analysis it is important to determine the 
entanglement Hamiltonian, which in this case has 
a single eigenvalue of zero. Extending this analy- 
sis to other single particle Hamiltonians, it can be 
shown for all reasonable single particle wave func- 
tions that the entanglement Hamiltonian will have 
an eigenvalue of zero if the chosen regions divide the 
probability density equally. 

If we model a bond with any single particle model 
in which probability density is split in half, the 
signature result is a zero eigenvalue of the entan- 
glement Hamiltonian. If one takes the associated 
Renyi entropy value, hi(2), as a normalization of 
the number of electrons being covalently shared be- 
tween two subsystems, then we can estimate the 
number of covalent bonds in the system with Nbond 
= S2/(2-ln(2)). Conventionally a bond is occupied 
by two electrons, which is why a factor of two is in- 
cluded in this definition. For C2 and N 2 at their 
respective equilibrium geometries this is suggestive 
of three covalent bonds. This is not the standard 
two bonds that one would expect from traditional 
bond order theory of C 2 [27-29]. We note that the 
bond order of C2 is still an open question and has 
been the subject of recent studics[30, 31]. As for the 
other dimcrs, the results for Li2 and F2 arc in agree- 
ment with conventional bond order values. However 
Be2 is well known to have unusual bonding charac- 
teristics, and its conventional bond order is zero[32]. 
This is problematic as Be2 is known to be a stable 
molecule and our result of a partial bond for Be2 
is what should be expected from physical consider- 
ations. 

It is important to further extend this analysis as 
the Renyi entropy only provides an average over 
some of the information in the density matrix pA- 
In particular not all contributions to the Renyi en- 
tropy are necessarily due to perfectly covalent bonds. 
More insight can be gained by considering the eigen- 
values of the entanglement Hamiltonian [33]. For the 
single determinant Hartree-Fock wave functions the 
eigenvalue spectrum can be calculated directly with 
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the correlation method. In this case the 2^ eigen- 
values of pa are mapped on to the N eigenvalues of 
the entanglement Hamiltonian. 

We proceed by observing that in the bonding 
region of C2 there is good agreement between 
the Hartree-Fock and multi-determinant Renyi en- 
tropies, as seen in figure 3. This indicates that the 
spectrum produced by the Hartree-Fock entangle- 
ment Hamiltonian may give an accurate represen- 
tation of the multi-determinant wave function spec- 
trum. The Hartree-Fock entanglement Hamiltonian 
for C2 and N2 contain several zero eigenvalues. For 
C2 there are 4 zero eigenvalues, and for N2 there 
are 6 zero eigenvalues, which is what one would ex- 
pect in traditional measurements of bond order for 
these molecules. The eigenvalue spectrum for C2 is 
plotted in figure 4. These two systems differ how- 
ever in terms of the remaining eigenvalues. In the 
case of N2 the rest of the eigenvalues are large and 
do not contribute significantly to the entanglement, 
whereas in C2 there are four more non-trivial values 
in the eigenvalue spectrum as shown in figure 4 by 
the two curves with diamond symbols. 

Our results suggest that there are two covalent 
bonds for C2 and three covalent bonds for N2. The 
other valence electrons in C2 appear to be strongly 
covalent as they make significant contributions to 
the Renyi entropy. It is likely these eigenvalues are 
a signature of the two extra bonds recently predicted 
in a study on bonding in C2PO]. In this study an 
energy was determined for each of the four bonds in 
C2 resulting in three energetically strong bonds and 
one weak bond. This is in contrast to our entan- 
glement results which show two perfectly covalent 
bonds and two partially covalent bonds. Differences 
between these results likely arise from our use of 
the entanglement Hamiltonian for the Hartree-Fock 
wave function. It is possible to determine the eigen- 
values of pa for multi-determinant wave functions 
through the calculation of many higher order Renyi 
entropies, S„[16]. However, this requires significant 
computational resources for C2 and developing more 
efficient techniques for this purpose is certainly an 
interesting direction for future research. 

In conclusion we report one of the first ab ini- 
tio calculations of Renyi entropy including Coulomb 
interactions while demonstrating the effect of inter- 
actions on proposed experimental measurements of 
entanglement. In general the fluctuation and corre- 
lation methods do not appear to accurately deter- 
mine the Renyi entropy. However for the systems 
under consideration here these methods bracket the 
Renyi entropy. 

Intuitively the use of entanglement appears to be 
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FIG. 4. Eigenvalues of the Hartree-Fock entanglement 
Hamiltonian calculated for the C2 dimer. There are 
twelve eigenvalues total, all non-zero eigenvalues are dou- 
bly degenerate and the zero eigenvalue is quadruply de- 
generate. Eigenvalues with large magnitude contribute 
less to the entanglement and the four eigenvalues with 
the largest magnitudes do not contribute to the entangle- 
ment significantly at any point along the binding curve. 



a natural way to characterize bonds and our re- 
sults suggest there are many possibilities for apply- 
ing these techniques to chemical systems. We have 
shown that entanglement can be used to help de- 
scribe covalent bonding in molecules and may be 
useful to describe other bonding situations. Other 
molecular systems of interest would be ionic and van 
der Waals bonded systems as well as transition metal 
dimers where the bonding properties can be difficult 
to characterize. 

In the case of more complicated molecules it will 
be important to consider the use of entanglement 
properties to generate partitions in space, similar 
to what is done with Bader analysis [34]. Finally 
we note that that there is large interest in using 
density functional theory for calculating entangle- 
ment entropy, where methods are already being 
developed [9]. It is clear that the results presented 
in this work as well as those from future QMC stud- 
ies using the techniques described here can serve as 
benchmarks for such developments. 
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II. SUPPLEMENTAL MATERIAL 

We have used several different techniques to calcu- 
late the entanglement properties in this work. Here 
we describe our implementation and usage of these 
algorithms. 

A. Particle Fluctuation Method 

The fluctuation method is calculated with mo- 
ments (and the related cumulants), of the particle 
number operator, which are straightforward to esti- 
mate in a VMC calculation. In non-interacting sys- 
tems the relationship between the cumulants and the 
S2 Renyi entropy is given as follows 

K+l 

S 2 = lim V a n (K)C n (2) 

K— >oo ' 

n=l 

with the cumulants, C n , and a set of constants, 
a n {K), that depend on the cutoff used in the ex- 
pansion, K. Similar expressions also exist for the 
rest of the Renyi entropies, and expressions for the 
coefficients a n (K) have been determined [16]. The 
S2 Renyi entropy converges quickly with the num- 
ber of terms kept in our results. We calculate the 
first six cumulant terms for use in this expansion 
in our VMC calculations. For testing the cumulant 
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cutoff wc compared against the swap operator for 
the Hartree-Fock calculations discussed in this work. 
The number of terms we included is more than suf- 
ficient to attain agreement between the two results 
within the error bars of their estimates. 



B. Correlation Method 

Calculating correlation functions of a quadratic 
Hamiltonian is another approach to calculating en- 
tanglement entropy that has been used significantly 
for calculations on lattice models. These calcula- 
tions are performed by mapping the Hamiltonian of 
the full system to the entanglement Hamiltonian, 
which is defined as a quadratic Hamiltonian that 
produces the density matrix of the subsystem. 



PA 



ne 



-He 



(3) 



The values k is a normalization that we use to fix 
Tv{pa) = 1- The correlation matrix is defined as 



(4) 



The entanglement Hamiltonian is then given by the 
correlation matrix and the identity matrix I[21], 



Hent = ln[ 



I-C, 



(5) 



The eigenvalues of the entanglement Hamiltonian 
can be used to calculate the eigenvalues of pa by 
equation 3. 

We were able to efficiently apply these methods to 
Hartree-Fock wave functions, for which this method 
is exact, and for multi-determinant wave functions, 
in which the correlation matrices only approximately 
determine the entanglement properties. The corre- 
lation matrices of interest for these wave functions 
can be calculated from the natural orbitals by pro- 
jection onto a real space grid that exists only in the 
subregion of interest. The straightforward way of 
doing this involves diagonalizing C'ij on an Ixl ma- 
trix, where I is the number of grid points in region 
A. The matrix C,j has a rank much smaller than 
the number of grid points, and a method that scales 
significantly better with the size of the grid can be 
performed by generating the following matrix from 
the natural orbitals m (rj ) and a grid rj in region A, 



Mi.. 



fl /2 Ui(rj). 



(6) 



The correlation matrix is formed by matrix multipli- 
cation MM T . Instead of diagonalizing this matrix, 



one can calculate the singular value decomposition 
of Mij. By squaring the resulting singular values 
one recovers the eigenvalues of the correlation ma- 
trix. This has the advantage of reducing the com- 
putational scaling with the number of grid points. 
We used a fine cubic grid of a minimum of I mil- 
lion points to minimize errors, and the sum of the 
eigenvalues of pa was normalized to one. 

The fi in the case of Hartree-Fock are equal to 1 
for occupied orbitals and for unoccupied orbitals. 
More generally they are equal to the eigenvalues of 
the natural orbitals for a given wave function. This 
technique can be used with any wave function based 
method, in which the natural orbitals can be calcu- 
lated. The natural orbitals can also be calculated 
within QMC, although this was not our approach in 
this work. 

It has been shown previously that only a sin- 
gle sum over the eigenstates of the entanglement 
Hamiltonian is needed to calculate the entanglement 
entropy [21]. A similar iterative formula over all the 
eigenstates can also be used to calculate the Rcnyi 
entropy S n . This is done by calculating a set of 
dummy variables w n and K n . We start by setting 
w a = and k — 1 and we calculate w i+ i and n i+ i 
iteratively, 



Wi+i = e ™ e, (l + Wi) + w t 



(7) 



After calculating WN tot and «JV tot j the Renyi entropy 
is given by 



WNtot + 1 



(8) 



The value n is the order of the Rcnyi entropy cal- 
culated and N tot is the number of eigenvalues. The 
number of non-trivial eigenvalues will at most be 
equal to the number of electrons in the case of the 
Hartree-Fock wave functions, or in the more general 
case it will be no larger than the number of natural 
orbitals used to calculate the correlation matrix. 

This brings up an ambiguity of applying the corre- 
lation method to a wave function that is not a single 
determinant. In principle the ground state correla- 
tion functions can be measured experimentally, and 
will agree with the correlation functions we calculate 
for high quality wave functions. Using such a corre- 
lation function yields a number of eigenvalues for the 
correlation matrix that exceeds the number of elec- 
trons in the initial system. In this work we applied 
equations 7 and 8 to the full set of non-trivial eigen- 
values that come from such a correlation matrix, but 
this is not the only approach for calculating Rcnyi 
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entropies when dealing with a correlation matrix for 
an interacting system. 



C. Swap operator method in quantum Monte 
Carlo 

Finally we extended a recent QMC method to cal- 
culate Renyi entropies exactly from a given wave 
function, to work on our molecular ab initio Hamil- 
tonians. A recent series of papers have shown that 
using a swap operator within lattice QMC calcula- 
tions yields a technique in which Renyi entropies can 
be calculated [12-1 5]. The operation of the swap op- 
erator is defined on a wave function that is written in 
terms of complete basis sets of regions A and B, \a) 
and |/3), as follows |*) = ^ZC at p\a)\0). The swap 
operator is defined in a space of the wave function 
that is a tensor product with itself as follows 

SWap A (^2 Cai,/3iK)|/3l» ® ( E A^.foK)!/^)) 

= E C «,ft E A, 2 ,fe(l"2>|/3i>®|ai>|/?2>). 



ot\,P\ a 2 ,P2 



(9) 



The expectation value of the swap operator is re- 
lated to pa by the following, 



2 = (* T ®V T \swap A \* T ®* T ) 



This operator can be sampled with pairs of walkers 
independently sampled from \^t\ 2 and the results 
can be used to calculate the Renyi entropy, S 2 = 
—ln(Tr(p a ) 2 ). In terms of the walker coordinates 
the expectation value of the swap operator is given 
by equation 11. From these definitions the VMC 
operator can be identified as 



0(a 1 ,a 2 ,/3 1 ,p 2 ) 



(10) 



ttr(fl(a2,ft))*r(fl(ai,/3 2 )) 
*T(fl(ai,/8i))*T0R(a 2 ,£2)) 

(12) 

In these equations we have used the notation 
R(a, f3), which are the coordinates of a walker. The 
swap operator is performed over two walker configu- 
rations, and thus the subscripts are used to identify 
whether the coordinates are from the first or second 
walker. All the coordinates in region A are associ- 
ated with a and all the coordinates in region B arc 
associated with (3. For our wave functions in which 
spin and particle number is held fixed, if the swap 
operator causes a change in spin or particle number, 
then the expectation values of the swap operator for 
that pair of configurations is zero. 
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(* T ® ^ T \swap A \^T ® *t) = y dxidx 2 • • • dx 2J v*T(-R(ai, /3i))*r(-R(a2, /3 2 ))*T(-R(a 2 , ^i))* T (-R(ai, &)) 

. / ^...^i^. W )Pi^.w>^g:;g;i^g;(u> 
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